home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / ctgsja.z / ctgsja
Encoding:
Text File  |  2002-10-03  |  11.4 KB  |  331 lines
  1.  
  2.  
  3.  
  4. CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))                                                          CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      CTGSJA - compute the generalized singular value decomposition (GSVD) of
  10.      two complex upper triangular (or trapezoidal) matrices A and B
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE CTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB, TOLA,
  14.                         TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
  15.                         NCYCLE, INFO )
  16.  
  17.          CHARACTER      JOBQ, JOBU, JOBV
  18.  
  19.          INTEGER        INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, NCYCLE, P
  20.  
  21.          REAL           TOLA, TOLB
  22.  
  23.          REAL           ALPHA( * ), BETA( * )
  24.  
  25.          COMPLEX        A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ), V(
  26.                         LDV, * ), WORK( * )
  27.  
  28. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  29.      These routines are part of the SCSL Scientific Library and can be loaded
  30.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  31.      directs the linker to use the multi-processor version of the library.
  32.  
  33.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  34.      4 bytes (32 bits). Another version of SCSL is available in which integers
  35.      are 8 bytes (64 bits).  This version allows the user access to larger
  36.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  37.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  38.      only one of the two versions; 4-byte integer and 8-byte integer library
  39.      calls cannot be mixed.
  40.  
  41. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  42.      CTGSJA computes the generalized singular value decomposition (GSVD) of
  43.      two complex upper triangular (or trapezoidal) matrices A and B. On entry,
  44.      it is assumed that matrices A and B have the following forms, which may
  45.      be obtained by the preprocessing subroutine CGGSVP from a general M-by-N
  46.      matrix A and P-by-N matrix B:
  47.  
  48.                   N-K-L  K    L
  49.         A =    K ( 0    A12  A13 ) if M-K-L >= 0;
  50.                L ( 0     0   A23 )
  51.            M-K-L ( 0     0    0  )
  52.  
  53.                 N-K-L  K    L
  54.         A =  K ( 0    A12  A13 ) if M-K-L < 0;
  55.            M-K ( 0     0   A23 )
  56.  
  57.                 N-K-L  K    L
  58.         B =  L ( 0     0   B13 )
  59.            P-L ( 0     0    0  )
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))                                                          CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))
  71.  
  72.  
  73.  
  74.      where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular upper
  75.      triangular; A23 is L-by-L upper triangular if M-K-L >= 0, otherwise A23
  76.      is (M-K)-by-L upper trapezoidal.
  77.  
  78.      On exit,
  79.  
  80.             U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R ),
  81.  
  82.      where U, V and Q are unitary matrices, Z' denotes the conjugate transpose
  83.      of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are
  84.      ``diagonal'' matrices, which are of the following structures:
  85.  
  86.      If M-K-L >= 0,
  87.  
  88.                          K  L
  89.             D1 =     K ( I  0 )
  90.                      L ( 0  C )
  91.                  M-K-L ( 0  0 )
  92.  
  93.                         K  L
  94.             D2 = L   ( 0  S )
  95.                  P-L ( 0  0 )
  96.  
  97.                     N-K-L  K    L
  98.        ( 0 R ) = K (  0   R11  R12 ) K
  99.                  L (  0    0   R22 ) L
  100.  
  101.      where
  102.  
  103.        C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
  104.        S = diag( BETA(K+1),  ... , BETA(K+L) ),
  105.        C**2 + S**2 = I.
  106.  
  107.        R is stored in A(1:K+L,N-K-L+1:N) on exit.
  108.  
  109.      If M-K-L < 0,
  110.  
  111.                     K M-K K+L-M
  112.          D1 =   K ( I  0    0   )
  113.               M-K ( 0  C    0   )
  114.  
  115.                       K M-K K+L-M
  116.          D2 =   M-K ( 0  S    0   )
  117.               K+L-M ( 0  0    I   )
  118.                 P-L ( 0  0    0   )
  119.  
  120.                     N-K-L  K   M-K  K+L-M
  121.  
  122.                M-K ( 0     0   R22  R23  )
  123.              K+L-M ( 0     0    0   R33  )
  124.  
  125.      where
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))                                                          CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      C = diag( ALPHA(K+1), ... , ALPHA(M) ),
  141.      S = diag( BETA(K+1),  ... , BETA(M) ),
  142.      C**2 + S**2 = I.
  143.  
  144.      R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
  145.          (  0  R22 R23 )
  146.      in B(M-K+1:L,N+M-K-L+1:N) on exit.
  147.  
  148.      The computation of the unitary transformation matrices U, V or Q is
  149.      optional.  These matrices may either be formed explicitly, or they may be
  150.      postmultiplied into input matrices U1, V1, or Q1.
  151.  
  152.  
  153. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  154.      JOBU    (input) CHARACTER*1
  155.              = 'U':  U must contain a unitary matrix U1 on entry, and the
  156.              product U1*U is returned; = 'I':  U is initialized to the unit
  157.              matrix, and the unitary matrix U is returned; = 'N':  U is not
  158.              computed.
  159.  
  160.      JOBV    (input) CHARACTER*1
  161.              = 'V':  V must contain a unitary matrix V1 on entry, and the
  162.              product V1*V is returned; = 'I':  V is initialized to the unit
  163.              matrix, and the unitary matrix V is returned; = 'N':  V is not
  164.              computed.
  165.  
  166.      JOBQ    (input) CHARACTER*1
  167.              = 'Q':  Q must contain a unitary matrix Q1 on entry, and the
  168.              product Q1*Q is returned; = 'I':  Q is initialized to the unit
  169.              matrix, and the unitary matrix Q is returned; = 'N':  Q is not
  170.              computed.
  171.  
  172.      M       (input) INTEGER
  173.              The number of rows of the matrix A.  M >= 0.
  174.  
  175.      P       (input) INTEGER
  176.              The number of rows of the matrix B.  P >= 0.
  177.  
  178.      N       (input) INTEGER
  179.              The number of columns of the matrices A and B.  N >= 0.
  180.  
  181.      K       (input) INTEGER
  182.              L       (input) INTEGER K and L specify the subblocks in the
  183.              input matrices A and B:
  184.              A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N) of A
  185.              and B, whose GSVD is going to be computed by CTGSJA.  See Further
  186.              details.
  187.  
  188.      A       (input/output) COMPLEX array, dimension (LDA,N)
  189.              On entry, the M-by-N matrix A.  On exit, A(N-K+1:N,1:MIN(K+L,M) )
  190.              contains the triangular matrix R or part of R.  See Purpose for
  191.              details.
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))                                                          CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))
  203.  
  204.  
  205.  
  206.      LDA     (input) INTEGER
  207.              The leading dimension of the array A. LDA >= max(1,M).
  208.  
  209.      B       (input/output) COMPLEX array, dimension (LDB,N)
  210.              On entry, the P-by-N matrix B.  On exit, if necessary, B(M-
  211.              K+1:L,N+M-K-L+1:N) contains a part of R.  See Purpose for
  212.              details.
  213.  
  214.      LDB     (input) INTEGER
  215.              The leading dimension of the array B. LDB >= max(1,P).
  216.  
  217.      TOLA    (input) REAL
  218.              TOLB    (input) REAL TOLA and TOLB are the convergence criteria
  219.              for the Jacobi- Kogbetliantz iteration procedure. Generally, they
  220.              are the same as used in the preprocessing step, say TOLA =
  221.              MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS.
  222.  
  223.      ALPHA   (output) REAL array, dimension (N)
  224.              BETA    (output) REAL array, dimension (N) On exit, ALPHA and
  225.              BETA contain the generalized singular value pairs of A and B;
  226.              ALPHA(1:K) = 1,
  227.              BETA(1:K)  = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = diag(C),
  228.              BETA(K+1:K+L)  = diag(S), or if M-K-L < 0, ALPHA(K+1:M)= C,
  229.              ALPHA(M+1:K+L)= 0
  230.              BETA(K+1:M) = S, BETA(M+1:K+L) = 1.  Furthermore, if K+L < N,
  231.              ALPHA(K+L+1:N) = 0
  232.              BETA(K+L+1:N)  = 0.
  233.  
  234.      U       (input/output) COMPLEX array, dimension (LDU,M)
  235.              On entry, if JOBU = 'U', U must contain a matrix U1 (usually the
  236.              unitary matrix returned by CGGSVP).  On exit, if JOBU = 'I', U
  237.              contains the unitary matrix U; if JOBU = 'U', U contains the
  238.              product U1*U.  If JOBU = 'N', U is not referenced.
  239.  
  240.      LDU     (input) INTEGER
  241.              The leading dimension of the array U. LDU >= max(1,M) if JOBU =
  242.              'U'; LDU >= 1 otherwise.
  243.  
  244.      V       (input/output) COMPLEX array, dimension (LDV,P)
  245.              On entry, if JOBV = 'V', V must contain a matrix V1 (usually the
  246.              unitary matrix returned by CGGSVP).  On exit, if JOBV = 'I', V
  247.              contains the unitary matrix V; if JOBV = 'V', V contains the
  248.              product V1*V.  If JOBV = 'N', V is not referenced.
  249.  
  250.      LDV     (input) INTEGER
  251.              The leading dimension of the array V. LDV >= max(1,P) if JOBV =
  252.              'V'; LDV >= 1 otherwise.
  253.  
  254.      Q       (input/output) COMPLEX array, dimension (LDQ,N)
  255.              On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually the
  256.              unitary matrix returned by CGGSVP).  On exit, if JOBQ = 'I', Q
  257.              contains the unitary matrix Q; if JOBQ = 'Q', Q contains the
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))                                                          CCCCTTTTGGGGSSSSJJJJAAAA((((3333SSSS))))
  269.  
  270.  
  271.  
  272.              product Q1*Q.  If JOBQ = 'N', Q is not referenced.
  273.  
  274.      LDQ     (input) INTEGER
  275.              The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ =
  276.              'Q'; LDQ >= 1 otherwise.
  277.  
  278.      WORK    (workspace) COMPLEX array, dimension (2*N)
  279.  
  280.      NCYCLE  (output) INTEGER
  281.              The number of cycles required for convergence.
  282.  
  283.      INFO    (output) INTEGER
  284.              = 0:  successful exit
  285.              < 0:  if INFO = -i, the i-th argument had an illegal value.
  286.              = 1:  the procedure does not converge after MAXIT cycles.
  287.  
  288. PPPPAAAARRRRAAAAMMMMEEEETTTTEEEERRRRSSSS
  289.      MAXIT   INTEGER
  290.              MAXIT specifies the total loops that the iterative procedure may
  291.              take. If after MAXIT cycles, the routine fails to converge, we
  292.              return INFO = 1.
  293.  
  294.              Further Details ===============
  295.  
  296.              CTGSJA essentially uses a variant of Kogbetliantz algorithm to
  297.              reduce min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and
  298.              L-by-L matrix B13 to the form:
  299.  
  300.              U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,
  301.  
  302.              where U1, V1 and Q1 are unitary matrix, and Z' is the conjugate
  303.              transpose of Z.  C1 and S1 are diagonal matrices satisfying
  304.  
  305.              C1**2 + S1**2 = I,
  306.  
  307.              and R1 is an L-by-L nonsingular upper triangular matrix.
  308.  
  309. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  310.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  311.  
  312.      This man page is available only online.
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.